Superfluidity of fermions with repulsive on-site interaction in an anisotropic optical 

lattice near a Feshbach resonance 
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We present a numerical study on ground state properties of a one-dimensional general Hubbard 
model (GHM) with particle assisted tunnelling rates and repulsive on-site interaction (positive-U), 
which describes fermionic atoms in an anisotropic optical lattice near a wide Feshbach resonance. 
For our calculation, We utilize the time evolving block decimation (TEBD) algorithm, which is an 
extension of the density matrix renormalization group and provides a well controlled method for 
one-dimensional systems. We show that the positive-U GHM, when hole doped from half-filling, 
exhibits a phase with coexistence of quasi-long-range superfiuid and charge-density-wave orders. 
This feature is different from the property of the conventional Hubbard model with positive-U, 
indicating the particle assisted tunnelling mechanism in GHM brings in qualitatively new physics. 
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The combination of Feshbach resonance and optical 
lattice techniques has opened up possibilities to investi- 
gate strongly interacting ultracold atoms under tunable 
configurations [1]. Ability to control such strongly inter- 
acting systems provides an unprecedented opportunity 
to explore interesting states of matter. Many interesting 
physics has been predicted for ultracold atom systems 
with fundamental Hubbard model hamiltonians. For ex- 
ample, with Bose-Hubbard model and its derivations, 
people have studied superfiuid to Mott-insulator transi- 
tion 0], existence of supersolid orderQ, etc, for ultracold 
bosons while with Fermi-Hubbard model, Luther-Emery 
Q and FFLO 0] phases are predicted to be observable for 
ultracold Fermions with attractive interaction. In partic- 
ular, it is well known that for repulsive (positive-?/) con- 
ventional (Fermi-) Hubbard model the susceptibility for 
superfiuid and charge density wave (CDW) orders are 
suppressed at low temperature and the leading quasi- 
long-range order is given by a spin density wave (SDW) 
at any filling fraction 

However, in this work, we show that coexistence of 
quasi-long-range superfiuid and CDW orders can be ob- 
served for fermionic atoms with repulsive on-site inter- 
action in an anisotropic optical lattice near a wide Fesh- 
bach resonance. The interactions in this strongly inter- 
acting system is described by a one-dimensional positive- 
U general Hubbard model (GHM) with particle assisted 
tunnelling rates The GHM is an effective one-band 
Hamiltonian that takes into account the multi-band pop- 
ulations and the off-site atom-molecule couplings in an 
optical lattice near a wide Feshbach resonance (see the 
detailed derivation in Ref. [7]). It is interesting to note 
that the GHM with similar particle assisted tunnelling 
also arises in different physical contexts, as proposed 
in Ref. [1]. In contrast with the case of conventional 
positive-U Hubbard model, we show that the superfiuid 
and CDW emerge as dominant quasi-long range orders 
over spin orders for the positive-U GHM when the system 
is significantly hole-doped below half-filling, although at 



or very close to half-filling, the dominant correlation in 
GHM is still anti-ferromagnetic. This feature indicates 
that the particle assisted tunnelling in GHM brings in 
qualitatively new physics. It makes the effective inter- 
action in GHM doping dependent, showing different be- 
haviors with a possible phase transition in between. We 
get our results through numerical calculation based on 
the time evolving block decimation (TEBD) algorithm 
0j fiot - which, as an extension of the density matrix 
renormalization group method [lH, is a well controlled 
approach to deal with one-dimensional systems. We com- 
pare our numerical results with some known exact results 
for the conventional Hubbard model, and the remarkably 
precise agreement shows that the calculation here can 
make quantitatively reliable predictions. 

As shown in Ref. Q, a generic Hamiltonian to de- 
scribe strongly interacting two-component fermions in an 
optical lattice (or superlattice) is given by the following 
general Hubbard model: 



(1) 



where riia- = al^aia, rii = riij -I- Uii, fi is the chemical 
potential, denotes the neighboring sites, and al^ is 
the creation operator to generate a fermion on the site 
i with the spin index a. The symbol a stands for (4, |) 
for a = (Tii)- The Sg and St terms in the Hamilto- 
nian represent particle assisted tunnelling, for which the 
inter-site tunnelling rate depends on whether there is an- 
other atom with opposite spin on these two sites. The 
particle assisted tunnelling comes from the multi-band 
population and the off-site atom-molecule coupling for 
this strongly interacting system 0]. For atoms near a 
wide Feshbach resonance with the average filling num- 
ber {rii) < 2, each lattice site could have four different 
states, either empty (with state |0)), or a spin | or J, 
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atom (ajg.|0)), or a dressed molecule (rf||0)) which is com- 
posed by two atoms with opposite spins. The two atoms 
in a dressed molecule can distribute over a number of 
lattice bands due to the strong on-site interaction, with 
the distribution coefficient fixed by solving the single- 
site problem. One then can mathematically map the 
dressed molecule state d- 10) to a double occupation state 
a||a|||0) by using the atomic operators a . After this 
mapping, the effective Hamiltonian is transformed to the 
form of Eq. (1). The GHM in Eq. (1) reduces to the 
conventional Hubbard model when the particle assisted 
tunnelling coefficients Sg and St approaching zero, as one 
moves far away from the Feshbach resonance. Near the 
resonance, Sg and St can be significant compared with the 
atomic tunnelling rate t due to the renormalization from 
the multi-band populations and the direct neighboring 
couphng Q. 

We consider in this work an anisotropic optical lat- 
tice for which the potential barriers along the x, y di- 
rections are tuned up to completely suppress tunnelling 
along those directions. The system becomes a set of in- 
dependent one-dimensional chains. We thus solve the 
GHM in one dimension through numerical analysis. For 
this purpose, first we transfer all the fermion operators 
to the hard core boson operators through the Jordan- 
Wigner transformation 1^ . In the one-dimensional case, 
we can get rid of the non-local sign factor, and after 
the transformation the hard core boson operators sat- 
isfy the same Hamiltonian as Eq. (1). On each site we 
then have two hard core boson modes which are equiv- 
alent to a spin-3/2 system with the local Hilbcrt space 
dimension d = A. We can therefore use the TEBD algo- 
rithm to solve this pseudo-spin system [§] . Similar to the 
density matrix renormalization group (DMRG) method 
the TEBD algorithm is based on the assumption 
that in the one-dimensional case the ground state = 
' ' ■ Si =1 Cii...i„ |«i • • • hi) of the Hamiltonian with 
short-range interactions can be written into the following 
matrix product form: 
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where f'^'I^^ denotes the matrix associated with site-s 
with the matrix dimension x- When x = 1, the assump- 
tion reduces to the mean-field approximation, and for a 
larger the matrix product state well approximates the 
ground state as it catches the right entanglement struc- 
ture for ID systems ^E^. To use the TEBD algo- 
rithm, we just start with an arbitrary matrix product 
state in the form of Eq. (2), and evolve this state with 
the Hamiltonian (1) in imaginary time through the prop- 
agator e~^*. The state converges to the ground state of 
the Hamiltonian pretty quickly. From the final ground 
state in the matrix product form, one can efficiently cal- 
culate the reduced density operator and various corre- 
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FIG. 1: The numerical result for the Hubbard model com- 
pared with some known exact result, (a) Ground state en- 
ergy as a function of U at half-filling (energy in the unit of 
t), where data points marked by solid dots are from the ex- 
act Bethe ansatz solution while those marked by pentagram 
are from our numerical program; (b) The relative error in the 
ground state energy; (c) Real-space spin correlation function 
at the filling fraction (rii) = 0.5 and U = 8t, compared with 
the asymptotic form in Eq. (3) (solid curve) with Kp = 0.62 
and A = 0.13. (d) Similar to (c), except that (rii) — 0.75, and 
the corresponding Kp = 0.60, A =0.19. 



lation functions. This calculation has a well controlled 
precision since at each time step to update the matrix 
product state, the Hilbert space truncation error can be 
suppressed by choosing an appropriate matrix dimension 
X 19]. In this calculation, we use the infinite lattice al- 
gorithm by assuming that the lattice is bipartite and the 
ground state has a translational symmetry for each sub- 
lattice . This allows us to directly calculate the system 
in the thermodynamic limit. 

To show that our calculation is capable of making re- 
liable predictions, we first test our results by compar- 
ing them with some known exact results of the Hub- 
bard model in certain cases. For the Hubbard model 
at half-filling {mi) = {rm) = 0.5, the ground state en- 
ergy per site is known to have the analytic expression 
E ^ -4. 4i+"xp(L"[}/2)] ill the thermodynamic hmit 
from the exact Bethe ansatz solution [13], where Jo and 
Ji are Bessel functions and we have chosen the tunnelling 
rate t as the energy unit. In Fig. 1 (a), we show our nu- 
merical results for the ground state energy of the Hamil- 
tonian (1) with Sg ^ St ^ 0, and one can see that it 
agrees very well with the exact energy of the Hubbard 
model in particular when U > t. The error is in general 
smaller than 10^'^ as shown in Fig. 1(b). In this and the 
following calculations, we choose the matrix dimension 
X = 40. We have tried larger x which gives better pre- 
cision, but we choose x = 40 to have a faster speed and 
its precision is enough for our purpose. 

We have also tested the final state from our calculation 
by comparing its correlation functions with some known 
results. It is difficult to get correlations analytically from 
the Bethe ansatz solution, but from the bosonization ap- 
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FIG. 2: The numerical results for the spin {Sr), the CDW (-Dr), and the pair (Pr) correlation functions for the GHM with 
different particle assisted tunnelling rates and at different filling fractions, where 5g = Q 5t = Q for solid curves, 5g — 3t St — —6t 
for dashed curves, Sg = 3t St = —3t for dotted curves, and Sg = 7t 5t = —lit for dash-dotted curves. 



proach to the one-dimensional Hubbard model, we know 
its correlation functions take certain asymptotic forms. 
For instance, one can look at the spin-spin correlation, 
defined as Sr = (s^ • Si+r), where the spin operator for 
the site i is given by Si = al^aa/3ai[3/2 with a and f3 —i, 
t and a standing for the Pauli matrices. The correlation 
Sr is independent of i because of the translational sym- 
metry. The Hubbard spin correlation function has the 
following asymptotic form 14 1 
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where Kp is the Luttinger parameter whose value has 
been determined before from the exact Bethe ansatz so- 
lution (isl. [l^. kp is the fcrmi momentum related to the 
filling number (ui) through kp = (rii) n/2, and A is a 
non-universal model dependent constant. In Fig.l (c) 
and (d), we compare our calculation results for Sr with 
this asymptotic form for filling number (n^) = 0.5 and 
0.75, and the agreement is again remarkable as long as r 
is not too small (the expression of Sr in Eq. (3) is not 
accurate for small r). 



With the confidence in numerics built from the above 
comparison, we now present our main calculation results 
for the repulsive GHM in Eq. (1) with J7 > 0. Apart from 
the spin correlation Sr defined before, we also calculate 
the charge-density-wave (CDW) correlation, defined as 
Dr = {rLirii^r) ~ {n-i) (rii-^-r) , and the pair (superfluid) 

correlation, defined as Pr = {o-ii CLiiO'l+ri^l+r\) ■ 
suits are shown in Fig. 2 for different filling fraction 
(rii) and for models with different particle assisted tun- 
nelling rates 6g and 6t. First at half filling with (n^) = 1, 
the correlation functions Sr, Dr, and Pr for the GHM 



with different dg and dt all look qualitatively similar to 
the corresponding results for the conventional Hubbard 
model, although with increase of the coefficient Sg the 
spin correlation reduces a bit while the CDW and super- 
fluid correlations increase slightly. Clearly, the dominant 
correlation in this case is in spin which suggests a quasi- 
long range anti-ferromagnetic order. In this and follow- 
ing calculations, we take U = 8t for all the cases, which 
corresponds to a significant on-site repulsion. 

Qualitatively different results show up when the sys- 
tem is doped with holes. At the filling fraction (n^) = 
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FIG. 3: The spin, the CDW, and the pair correlation functions 
in momentum space for the GHM with &g = J,t it = —?,t. The 
solid, dashed, and dash-dotted curves correspond to filling 
factor (jii) = 1, 0.75, and 0.5, respectively. For the calculation 
of the Fourier transformation, we have used the real space 
correlation functions for A*" = 100 sites. 

0.75, although for the Hubbard model the spin correla- 
tion is still the dominant one (the spin density wave or- 
der has been pinned to the corresponding 2kp = 37r/4), 
for the GHM with a noticeable the superfluid and the 
CDW emerge as the leading quasi-long-range orders, and 
their correlations increase significantly and decay much 
slower in space compared with the spin correlation when 
8g grows. These features become more evident when we 
further increase the doping. For instance, at the right col- 
umn of Fig. (2), we show the correlations for the filling 
fraction (n^) — 0.5. The qualitative behavior is similar 
to the case with (n^) = 0.75, but the CDW and super- 
fluid correlations for the GHM get significantly larger at 
long distance, and the contrast with the Hubbard model 
becomes sharper. One also note that for all these calcu- 
lations, change of the coefficient 8t in the GHM makes 
little difference to the result. This is understandable as a 
significant positive U suppresses the possibility of double 
occupation in the lattice, and the 8t term in the GHM 
has no effect without double occupation. The lig term in 
the GHM, however, is critically important, which favors 
superfluidity in general and brings in the qualitatively 
different features mentioned above. 

To show the spatial structure of these quasi-long range 
(QLR) orders, we plot in Fig. 3 the spin, the CDW, and 
superfluid correlations in the momentum space for the 
GHM with 8g = 3t at different filling fractions (rii). The 
momentum space correlations are defined by the Fourier 
transform Xk = J2r=o -^r cos{kr), where X stands 

for the correlations 5, D, or P. From these momentum 
space curves, one can clearly see that this GHM at half 
filling has a QLR anti-ferromagnetic order (characterized 
by the peak at /c = tt) , and away from half filling a QLR 
superfluid order (peak at fc = 0) and CDW order (peaks 



at fc = 2kF and 2tt — 2kp, where 2kF = 37r/4 {tt/2) for 
the fiUing fractions (n^) = 0.75 (0.5), respectively). The 
peaks in Fig. 3 have finite widths because these orders in 
ID are only quasi-long-range with algebraic decay. Note 
that if we turn on small tunnelling interaction between 
different ID chains, a leading QLR order, such as super- 
fluid order, could be stabilized to a true long range order 
. The GHM thus provides an example of a microscopic 
Hamiltonian that with hole doping from half filling, an 
anti-ferromagnetic phase could be transferred to a su- 
perfluid phase (or a CDW phase in some case depending 
on which order becomes more dominant with the inter- 
chain coupling). The correlations that characterize these 
QLR orders can be detected for the cold atomic gas, for 
instance, through the method described in Ref. fig . 

In summary, we have investigated the ground state 
properties of the general Hubbard model with repul- 
sive on-site interaction in one dimension through well 
controlled numerical analysis. For the system with sig- 
niflcant particle assisted tunneling rates 6g and St, we 
have found coexistence of quasi-long range superfluid 
and charge-density-wave orders when the system is hole- 
doped from half filling. This feature is in sharp con- 
trast with convention Hubbard model, in which case 
for positive-?/ the charge and superfluid orders are al- 
ways suppressed regardless of the flUing fraction. With a 
combination of the Bosonlization approach and the nu- 
merical method here, it may be possible to determine 
the compete phase diagram for the GHM. The model 
here describes strongly interacting fermionic atoms in an 
anisotropic optical lattice. The possibility of a transition 
from an anti-ferromagnetic phase to a superfluid phase 
for the GHM with hole doping may also have interesting 
indications for other areas. 
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